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Abstract 

A new class of high-order monotonicity-preserving schemes for the numerical solution of conserva¬ 
tion laws is presented. The interface value in these schemes is obtained by limiting a higher-order 
polynomial reconstruction. The limiting is designed to preserve accuracy near extrema, and to work 
well with Runge-Kutta time stepping. Computational efficiency is enhanced by a simple test that 
determines whether the limiting procedure is needed. For linear advection in one dimension, these 
schemes are shown to be monotonicity-preserving and uniformly high-order accurate. Numerical 
experiments for advection as well as the Euler equations also confirm their high accuracy, good 
shock resolution, and computational efficiency. 
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Introduction 


We consider higher-order schemes (at least third-order) for the numerical solution of the Euler 
equations. Typical solutions to these equations have smooth structures interspersed with disconti¬ 
nuities. The challenge is to develop schemes that are highly accurate in smooth regions and have 
sharp nonoscillatory transitions at discontinuities. 

Achieving this dual objective remains a daunting task. Among the first attempts, Colella and 
Woodward [2] introduced a piecewise parabolic method (PPM), which employs a four-point centered 
stencil to define the interface value; this value is then limited to control oscillations. The centered 
stencil, however, results in a scheme with a large dispersion error, and the limiting procedure causes 
accuracy to degenerate to first-order near extrema. 

The essentially nonoscillatory (ENO) schemes of Harten et al. [3] were developed via a different 
line of thought. In these schemes, an adaptive stencil is used to select the “smoothest” data, 
thereby avoiding interpolations across discontinuities. While an adaptive stencil does avoid spurious 
oscillations near discontinuities, it does not make use of all the available data. The weighted-ENO 
(WENO) schemes by Liu et al. [11] and Jiang and Shu [9] make better use of the available data 
by defining the interface value as a weighted average of the interface values from all stencils. The 
weights are designed so that in smooth regions the scheme nearly recovers a very accurate interface 
value using all stencils but, near discontinuities, it recovers the value from the smoothest stencil. 
The WENO schemes, however, are still diffusive: they smear discontinuities nearly as much as the 
ENO schemes. 

In this paper, we follow the limiting approach. The interface value is defined by a five-point 
stencil. As a result, the leading error is dissipative, and the dispersive error is considerably smaller 
than that of the four-point stencil. This interface value also combines well with the three-stage 
Runge-Kutta time stepping. Similar to PPM, oscillations an' controlled by a limiting procedure. 
The key differences, however, are that our limiting is designed to preserve accuracy near extrema 
and to work well with Runge-Kutta time stepping. The resulting scheme is accurate in smooth 
regions, resolves discontinuities with high resolution, and is also efficient. 

Note that a piecewise linear scheme of this type was presented by Huynh [(>]. Extensions to 
piecewise parabolic schemes were presented by Suresh [18] and Huynh [7]. The present scheme 
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incorporates these ideas within a Runge-Kutta time integration framework. 

In §1, the spatial discretization and the Runge-Kutta time integration are reviewed. Section 2 
describes the reconstruction procedure, which is the key feature of our scheme. Extensions of this 
scheme to systems of equations and multi-dimensions are dealt with in §3. Numerical experiments 
appear in §4. Finally, conclusions are presented in §5. 

1. Discretization 

For simplicity, we describe the methods for the advection equation with constant speed a, 

u t + au x — 0, (1.1 a) 

u(x, 0) - uq ( x ) (1.1 b) 

where / is time, x is distance, and Uo(x) is the initial condition. For the moment, the solution is 
assumed to be periodic in x so that no boundary conditions are needed. 

Let xj be the cell center of a uniform mesh, x J+[ / 2 the interface between the j-tli and j + 1-th 
cells, and h the cell width. Denote by Uj(t) the cell average of u at time /, 

Uj(t) = - / u(x,t)dx. (1.2) 

“ x J — \/2 

Integrating (1.1a) over the cell [x J _ 1 / 2 , ^+ 1 / 2 ] yields 

~~dt~ Ti u { x 3 + x ! 2 ' ^ — u ^ x J* 1 / 2 ’ ^ ~ (!•’*) 

At time t n = nr where r is the time step, assume that we know v 71 which approximates iij(t n ). 
We wish to calculate For simplicity of notation, we omit the superscript n when there is no 

confusion, e.g., Vj denotes v 77 . 

An approximation to u(xj+i / 2 ,£ n ) is called the interface value and is denoted by v J+l / 2 - The 
calculation of the interface value from the known cell averages is accomplished in two steps. In 
the first or reconstruction step, nonoscillatory approximations of / 2 ,/ n ) to the left and right 

sides of the interface Xj+i/ 2 denoted by t^ +1 y 2 and 1 } ^. x / 2 are constructed. This step determines 
the scheme's order of accuracy and is the main concern of this paper. In the next or upwind step, 
the interface value is determined bv the wind direction: If a > 0, Vj+ 1/2 - r? f+i/2’ ot-herwiso, 
Vj+xj 2 = v f+\/ 2 ' Thus for advection, we need only one of the two values and for the 
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Euler equations, however, we will need both, and we employ well-known methods for the upwind 
step. 

Equation (1.3) can be integrated by a standard Runge-Kutta method. Here we use the three- 
stage scheme of Shu and Osher [15]. With v representing denote by L(v) the spatial operator 

L(v)j = -(Vj+i /2 - v j-i/ 2 )- (1-4) 

Then this scheme is given by 

~ V n 

w( l ) = iy(°) -f aL(w^) 

- | w (°) + + aL{w^)) (1.5) 

-f |(u;( 2 ) + ctL{ w^)) 

jg* “f* 1 — yj (3) 

where a = arjh is the CFL number. 

Observe that Runge-Kutta schemes like (1.5) are made up of repeated applications of a single 
stage scheme given by w {k) + <r/„(«><*>). A- = 0. 1, and 2. Moreover, each stage is an explicit Euler 
scheme, e.g., 

w j = V J - <r(Vj+l/2 - Vj- 1/2 )• ( 1 .<>) 

Therefore, we first design a monotonicity-preserving scheme for (1.6) and then extend it to the full 
scheme (1.5). 


2. Reconstruction 

Without loss of generality, we discuss the reconstruction only for r^ +1/2 , i.e., we assume a > 0. 
The reconstruction is carried out in two steps. In the first step an accurate and stable formula is 
used to compute the interface value which is called the original value. In the second step, this value 
is then modified or limited appropriately to achieve a monotonicitv-preserving scheme. 

A straightforward choice for vf + i /2 usin § tllc f ' 1VG co11 averages i>j_ 2 . v J+2 (the same stencil 

as the third-order ENO scheme) is 

? >H/2 ~ - 13vj-\ + 4iVj + 27e J+1 l]Vj+ 2 )/ 60. (2.1) 

Oth er choices include a low phase error fourth-order formula [8] 

vf+i /2 = (9«>-2 - , + 191 Vj + 104» j+1 - 1 lvj+ 2 )/240, (2.2) 
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or a fifth-order accurate implicit formula given by 

(3i>jl,/2 + fit’j+1/2 + v j+3/2)/ 10 = ( W J-1 + 19i ’j + 10u J+ i)/30. (2.3) 

The implicit formula has the advantage of low dispersive and dissipative errors; its disadvantage is 
that the tridiagonal matrix inversion costs more. 


2.1 Monotonicity constraint 

We derive constraints for the interface value so that monotonicity is preserved by (1.6). First, 
we need a few definitions. Let the median of three numbers be the number that lies between the 
other two. Let minmod (x, y) be the median of x , y , and 0. Equivalently, 

minmod (x, y) = \ [sgn(.r) + sgn(j/)] min (|x|, | y \). (2.4) 

Conversely, the median function can be expressed in terms of minmod, 

median (x, y,z) — x + minmod (y — x,z — x). (2.5) 

The minmod function can be extended to any number of arguments. For k arguments, minmod ( z x .. . 
returns the smallest argument if all arguments are positive, the largest if all are negative, and zero 
otherwise. This function can be coded as 


minmod (z \,.. ., z^) - s min(|zi|. \zk \) 


(2.6a) 


where 


5 = i(sgn(zi) + sgn(z 2 )) ^(sgn(~i) + sgn(- 3 ))...i(sgn(c 1 ) + sgn(^.)) . (2.66) 

Also denote by I [z 1? ..., Zk] the interval [min(Ti,..., z*), max(:i,..., z^)]. 

We can now derive a simple condition that preserves monotonicity. At interface j - 1/2, suppose 
the value lies between Vj-\ and v 3 : 


V j- 1/2 ^ / [uj_i, Xj]. 


(2.7) 


Next, for the interface j + 1/2, denote 


V ( ' L = Vj + a(vj -Vj-i) 


( 2 . 8 ) 
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X j~ 1 X j X j+ 1 X 

Figure 2.1: Monotonicity-preserving constraint (2.7) and (2.9). 


where UL stands for upper limit, and a > 2 (more on a momentarily). Suppose the value v 1 ' , 

11 J + 1/2 

lies between Vj and v UL , 

v j+ 1/2 € I [vj, v UL }. (2.9) 

Then, after one stage via (1.6), the solution ie] 1) lies between r_,_i and v y provided that, the time 
step satisfies the condition 

<7 <1/(1 +a). (2.10) 

Indeed, for increasing data, (2.7) and (2.9) imply that the steepest slope v UL - v J _ l satisfies 
l ' L L ~ v j-i — ( Q + 1)( l 'j ~ v j- 1 ); therefore, (2.10) implies v J _ i < tr* 1 * < Vj. See Fig. 2.1. 

Note that for parabolic reconstruction schemes, a is typically 2. For Runge-Kutta time stepping, 
we find that a = 4 works well, while o = 2 tends to cause stair-casing. With « = 4, (2.10) leads to 
a GFL number o < 0.2 but, in practice, <7 = 0.4 still yields nonoscillatory results. 

Next, assume that (2.7) and (2.9) hold for all j. Expression (2.7) with index j replaced by j + 1 
takes the form 


V J + 1/2 € I[vj,Vj + ]]. 

The above and (2.9) result in the condition that v L . , lies 

] + 1/i 

/ {v 3 , v j+i] and v UL ], One end of this intersection is Vj 
and v UL , and is denoted by v MP . Using vj as the pivot, v MP 
function, 


( 2 . 11 ) 

in the intersection of the two intervals 
. The other is the median of vj. i^ +1 , 
can be expressed by using the minmod 


All 


— v 


+ minmod [?>j +1 - vj , o( 


- v 


j-i 


>]■ 


( 2 . 12 ) 
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Thus, (2.9) and (2.11) imply 

*4+1/2 ^ I[vj,v MP ], (2.13) 

The simplest way to satisfy this constraint is to replace the original ^ L +1 / 2 Kv ^ ie median of v F +l j 2 , 
Vj, and v MP : 

v j+ 1/2 median (vf +l/2 , Vj, v MP ). (2.14) 

Expression (2.14) preserves monotonicity in the following sense: under the CFL restriction 
(2.10) if the data { v 3 } are monotone, then after one stage, {u^} are also monotone. This fact 
follows because lies between v 3 _i and v 3 for all j. 

The monotonicity-preserving property extends easily to the full scheme (1.5). Indeed, given 
monotone data {^}, we have just shown that are monotone provided the interface values 

are given by (2.14) and the above CFL restriction is satisfied. Since {u^ 1 ^} are monotone, the 
quantities {(C 1 * _|_ aL(vM^)} are also monotone because they result from applying a single stage 
scheme to Next, for each j, wj 2 ' is a combination of vj and (u^ 1 ^ + a{L( w^ l ^)) J with positive 

weights independent of j. Therefore, {u^ 2 '} are monotone as well. Repeating this argument, it 
follows that {} are also monotone. Thus, if the data { vj } are monotone, and the interface 
values are obtained by (2.14), then the cell averages at the next time level {n" +1 } are monotone 
under the CFL restriction (2.10). 

The drawback of (2.14) is that near an extremum, it causes accuracy to degenerate to first- 
order. Figures 2.2(a) and 2.2(b) show the loss of accuracy caused by the constraints (2.11) and 
(2.9), respectively. Note that the data are on a parabola. 


2.2 Accuracy-preserving constraint 

To avoid the loss of accuracy, we enlarge the intervals in (2.11) and (2.9) in such a way that 
these intervals remain the same for monotone data but, near an extremum, these intervals are 
larger, and both contain the original v ^ij 2 - 

First, the interval in (2.11) is enlarged by adjoining the value v MD defined below (MD stands 
for median). At interface j + 1/2, let v FL and v FR be the values extrapolated linearly from the left 
and right, respectively. 


- Vj + \( Vj - Vj-i ), v f R - v 3 + \ + |( Vj + l “ Vj + 2 ). 


( 2 . 15 ) 


7 




With 


Figure 2.2: Loss of accuracy near extrema: (a ) by (2.11) and (b) by (2.9). 

e — 2 ( v j v 3 +* )• 


where AV stands for average, set 


(2.16) 


v 


Ml) 


= median (v A v , v 1 ' v FR ). 


(2.17) 


Constraint (2.11) is relaxed to 


,! t+i /2 e / [rj. v j+u v MD ] 


(2.18) 


One can verify that if the four pieces of data {, Vj, m+j, v J+2 } are monotone, then, at the 
interface j + 1/2, v MD lies between Vj and Vj +1 , and the above constraint reduces to (2.11). Near 
an extremum as in the case of Fig. 2.3 (a), however, v AfD lies outside I [vj, e J + I ] and provides room 
so that the interval I[vj, Vj +contains the original v^ x j 2 - 

The argument (2.15)-(2.18) conv r eys the idea, lor the purpose of coding, it is more efficient if 
we employ the second differences. Set 


d 3 - Vj _i + Uj + 1 - 2 vj. 


(2.19) 


and 


< l j + i/2 = tninmod </ J+1 ), 

where MM stands for minmod. Then, since v FL = v AV - \d, and similarly for v FR , it follows that 


2 J 


_ v AV 1 . . 

' - ~ 2 (l j + l/2- 
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(a) 


V 



x j+ 1/2 


(b) 


Figure 2.3: Enlarged monotonicity interval (a) by (2.18) and (b) by (2.23). 


Next, the interval in (2.9) is enlarged by adjoining the value v LC defined below (LC stands for 
large curvature). Consider the parabola p determined by the cell averages vj, and the second 

difference d (a quantity similar to d 3 ). A straightforward calculation gives the value at .r J + 1 / 2 < 

p(*j+l/2) = V J + \{ V J - V J -1 ) + l d - 

The parabola with d — 2gives 

v LC = Vj + Vj - Vj-1 ) + (2.22) 


Constraint (2.9) is relaxed to 


d+t / 2 et[v^v u \ v [ ' c ] 


(2.23) 


Using the fact that o > 2, one can verify that if the four pieces of data {rj_ 2 , i?j -\, vj , } 

are monotone, then, at the interface j + 1/2, v LC lies between Vj and v i L . The above constraint 
therefore reduces to (2.9). Near an extremum as in the case of Fig. 2.3 (b), however, v LC lies outside 
I[vj,v^L] and provides room so that the interval I [t; 2 , v^'] contains the original ■ 

In practice, we reduce the amount of room so that near a non-monotone discontinuity (such 
as a sawtooth profile), constraints (2.18) and (2.23) reduce to (2.11) and (2.9) respectively. This 
reduction can be accomplished by replacing dj+Y /2 by 


dj+\/ 2 — minmod (4d 3 — dj + i,4dj+i - dj,dj,d 3 + 1 ). 


(2.24) 


To clarify the role of d^ 4 , assume that d 3 and d 1+ ] are of the same sign. Then if d 3 /d 3 +\ < 1/1 
or f/j/f/j + i > 4, the above minmod of four arguments returns 0. Loosely put, when the second 
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differences change substantially, then d^ + \ /2 = 0, and the extended intervals reduce to the simple 
ones in (2.9) and (2.11). We thus replace expressions (2.21) and (2.22) for each interface j + 1/2 
by 

v MD = m' ,V ' - 
„LC .. , 1 / . 


la ted by 


v 


mm 


.max 


AID 

AID 


d M4 

2 + 1/2’ 

(2.25) 


(2.26) 

v J+ i,v MD ] and / [vj, v UL , v LC ] can 

be calcu- 

), min(t>j, v 1 ' 1 ', v LC )\, 

(2.27«) 

), max! Vj, v 1 lj , v lx ) ]. 

{2.27 h) 


Finally, to ensure that vj) +lf2 lies in [u min ,u max ], we replace vf +l/2 by the median of vf +l/2 , 
v min , and n max : 

v j+i /2 ~ median (n| +1/2 , n mm , t- max ). (2.28) 

The above limiting procedure preserves monotonicity and accuracy. In addition, for most cells 
in smooth regions, the original Vj +( ^ 2 satisfies constraint (2.13) a priori. In this case, the limiting 
procedure (2.24)-(2.28) does not alter the original vj' +l / 2 - As a result, we can use (2.13) to detect 
such cells and bypass the limiting procedure altogether. The condition that the original , 

° J 4-1 /2 

lies in the interval I[ Vj ,v MP ] is equivalent to (v]' +l/2 - Vj)(i >]' +l/2 - v MP ) < 0. In practice, this 
condition is coded with a tolerance value of c = 10~ 10 : 


( V j+l/2 - Vj)(Vj + 1/2 - V MP ) < (. 
We summarize the computation of the interface value below. 


(2.29) 


Algorithm for the interface value. Suppose the cell averages {tq} are given, and a > 0. For 
each interface j + 1/2, calculate the original value vf +1/1 from (2.1), and v MP from (2.12). If (2.29) 
holds, then i’j+ 1/2 = ai] d wc move on to the next interface. Otherwise, calculate dj-i ■ d,, 

dj +1 from (2.19), dj^j 2 and d'/j / 2 from (2.24), ipfrom (2.8), v A 1 from (2.16), v MD and v LC 
from (2.25) and (2.26), and v mm , e max from (2.27). Finally, calculate r| +I/2 from (2.28), and the 
interface value is then v J+l / 2 = 
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2.3 Remarks 

Note that constraint (2.13) on the interface value is a sufficient condition for monotone data 
to remain monotone under Runge-Kutta time stepping. It may be viewed as the analogue of Van 
Leer’s constraint [19] which provides the same type of condition for monotonicity under exact time 
evolution. Also note that the geometric framework, the use of the median function, and v MD were 
introduced by Huynh in [5]. 

For advection with a < 0, the interface value ^ X/ / 2 1S obtained by reflecting the above expres¬ 
sions about To be specific, the reconstruction algorithm is the same with {vj_ 2 , Uj, Vj+i, 

Vj +2 } replaced by {v J+ 3 , Vj+ 2 ,Vj+i, Vj,Vj -\} respectively. Next, the stencil for computing both 
v^i /2 an ^ v j + 1/2 cons i sts °f s * x points {vj_ 2 , ...Vj+ 3 }. Therefore, we could define both 
and v*_ l/2 by the quintic fit of all six cell averages without enlarging the stencil (in the case of 
Euler equations). The corresponding limited scheme, however, is prone to stair-casing. 

Higher-order schemes can be derived using larger stencils. With the same stencil as the /nth- 
order ENO scheme, a (2 m - l)th-order scheme can be obtained. For example, for m = 1. we have 


the seven-point formula 


>j + 1 / 2 - (-3^_ 3 + 25i?j- 2 - + 319vj 


+214i? J _i - 38i?j +2 + 4 u j+3 )/420, 


(2.30a) 


and, for m = 5, the nine-point formula 


f+1/2 = (4 vj- 4 - 41t7j_ 3 + mvj-2 - 6411V-1 + 18791-, 


(2.306) 


+ 1375v,+i - 305?;j+ 2 + 55 c j+ 3 - 5i>j+4 )/2520. 

The same limiting can be employed for these original interface values. The resulting schemes achieve 
high spatial accuracy but remain third-order in time. For the fourth- and fifth-order Runge-Kutta 
methods, in order to preserve monotonicity we need the calculation of the time-reversed operator 
L [15], which is beyond the scope of this paper. 

The above reconstruction depends continuously on the data in the sense that a small change in 
the data causes a small change in the interface value. This property is shared by WENO (but not 


by ENO) reconstruction. 


3. Extensions 


In this section, we describe the extensions of the above schemes to the fuller equations. While 
these extensions are standard, the monotonicity-preserving property may not hold because the equa- 
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tions are nonlinear. Nevertheless, the numerical solutions obtained below are generally nonoscilla- 
tory. 

3.1 Euler system in one dimension 

The Euler equations of gas dynamics for a polytropic gas can be written as 

u f + f(u) x = 0 (3.1) 

where 

u = (p,pu, E) T , 

f(u) = uu + ( 0 ,p, up) T , ( 3 - 2 ) 

P =(7~1 )(E-\pu 2 ). 

Here, T represents the transpose; />, u, p, and E are tlie density, velocity, pressure, and total energy 
respectively; and 7 = 1.4, is the ratio of specific heats. The speed of sound c is given by ( 7 p/p) 1 ? 2 . 

The eigenvalues of the Jacobian matrix A(u) = df/du are u - c. u and a + c. The matrices of 
left and right eigenvectors of A are needed in the reconstruction. These are given bv ([ 3 ]) 

( ^ 2/2 T w/2c —b\u/2 — 1/2 c b\j2 \ 

L = 1 b’2 biU -61 (3.3) 

^ 6 2 /2 - u/‘2c — bjit/2 + l/2r b { /2 / 

and 

1 1 1 I s ! 

R= u-c u u + c ( 3 ..J) 

^ II - uc \u 2 H + uc 

where bi = (7 - l)/c 2 and b 2 = u 2 b x / 2 , H = c 2 /(7 - 1 ) + \u 2 . 

Integrating (3.1) over the cell [.r J _ 1/2 , .^+ 1 / 2 ] yields 

du, 1 r 1 

~dT + h [ f(u(;i V+ 1 / 2 ,^))-f(u(^-i/ 2 . 0 )J = 0 ^ ( 3 . 5 ) 

where u j(t) are the cell averages. The first step in calculating f J+1/2 ~ f(u(.T J+1 / 2 , t n )) is to 
reconstruct u on both sides of the interface x J+l / 2 . 

It is well known that the reconstruction is best carried out in local characteristic variables [3], If 
{vj} are the approximations to the cell averages at time level n, these local characteristic variables 
for the cell [£j_i/ 2 i Ti+i/s) are given by 

w A . = L (Vj)v j+k , for k = -2,2. (3.0) 
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The scalar reconstruction algorithm is now applied to obtain point values w^ 2 and wf /2 at x JJrl j 2 
and xj_i/ 2 , respectively. The corresponding conservative variables v^ +1 ^ 2 and v^_^ 2 are obtained 
by the inverse of (3.6) 

v f + l/2 = R ( V j) W 1/2’ v f-l/2 = R ( V j) W -1/2- ( 3 - 7 ) 

At each interface j + 1/2, the two values and v^ +1/2 are used to calculate f J+1 / 2 v * a Roe’s 

flux-difference splitting [12]. This splitting is implemented here with Huynh’s entropy fix [6]. 

Equation(3.5) is then integrated by the Runge-Kutta scheme (1.5). The time step is given in 
terms of the CFL number a by 

Af = - °—4 -• (3.8) 

Note that the extension described above is standard, but it does not take advantage of the fact 
that our reconstruction algorithm leaves the left and right interface values unchanged in smooth 
regions away from extrema. In these regions, the reconstruction applied to the local characteristic 
variables yields a result identical to formula (2.1) applied to the conserved variables {v ; }. Thus, the 
expense of characteristic decomposition may be avoided for such regions if they could be detected 
in a simple manner as in [6]. However, we do not pursue this approach here. 


3.2 Euler system in two dimensions 

An immediate extension of the above scheme to multi-dimensions can be accomplished in the 
same manner as the finite difference ENO schemes of Shu and Osher [15] [16]. The idea is to 
avoid calculating the mixed derivatives of the reconstruction from cell averages by applying the 
reconstruction directly on point values of the fluxes. The reconstruction then reduces to two one- 
dimensional reconstructions along the coordinate lines. The same Runge-Kutta scheme (1.5) is used 
to integrate the equations in time. Here we have chosen the Lax-Friedrichs version (ENO-LLF) in 
our numerical experiments. The 2D extension of the schemes derived here are obtained by substi¬ 
tuting our algorithms for reconstruction in place of the scalar one-dimensional ENO reconstruction 
therein. Coding aspects of these schemes can be found in [17]. 

4. Numerical experiments 


For simplicity, we present numerical results only for the scheme combining the quartic fit (2.1) 


and the accuracy-preserving constraint (2.2K). 


We refer to this scheme as the MP5 scheme (MP 



for monotonicity preserving). Some comparisons with EN03 and WEN05 schemes are provided. 
The three schemes MP5, WEN05 and EN03 have the same stencil, and the first two are spatially 
fifth-order accurate. Listings of the these three reconstruction procedures in Fortran are given in 
Appendix A. Also note that we employ only uniform meshes and, unless otherwise stated, the CEL 
number is 0.4. 

All computations are carried out on a 100 MHz R4000 SGI Indigo Workstation, with c = 10~ 10 
and a = 4. In all cases, the compiler options -r8 -03 were used. We have observed that computing 
times vary widely depending on the hardware and compiler options used. Therefore, computing 
times are to be viewed only as an approximate measure of the efficiency of the various schemes. 
The computing time for the scheme with constant reconstruction and the three-stage Runge-Kutta 
time stepping is also provided as a reference. Since the reconstruction is trivial for this scheme, 
this computing time reflects the cost of all other calculations except reconstruction. 


4.1 Advection of a smooth function 

We solve (1.1) with n(x,0) = sin(7rx) 4 with periodic boundaries. We are particularly interested 
in the behavior of the errors of the cell averages under mesh refinement. Since the function is 
smooth, the most accurate and efficient scheme with the given stencil of five cell averages is the 
unlimited scheme (2.1). We compare the results of the WEN05 and MP5 schemes to this unlimited 
scheme for a = 0.05 in Table 2 (a) and rr = 0.1 in Table 2 (b). The results from EXO 3 are loss 
accurate and are not shown. 

Note that the errors obtained by the unlimited scheme and those by the MP5 scheme are 
essentially identical. This confirms that the limiting procedure leaves the quartic fit unchanged 
at smooth extrema. At low CFL numbers, the MP5 scheme approaches the theoretical order of 
accuracy of five as can be seen in fable ‘2(a). In both cases, the MP5 scheme compares favorably 
with the WEN05 scheme in both accuracy and efficiency. 



Table - 2: Advection of sin(7rx) 4 by several schemes, t = 2, Ax = 2/iV, CPU time quoted is for 


all grids. 


Table 2(a): At/Ax = 0.05, CPU time for constant reconstruction: 4.30 sec. 


Scheme 

N 

error 

L t oc> order 

L\ error 

L\ order 

CPU time - sec. 


16 

2.39(-l) 

- 

1.07(1) 

- 



32 

3.45(-2) 

2.79 

1.73(-2) 

2.62 


WEN05 

64 

3.51(-3) 

3.29 

1.75(-3) 

3.31 

22.33 


128 

3.44(-4) 

3.35 

8.88(-5) 

4.30 



256 

1.15(-5) 

4.90 

2.54(-6) 

5.13 



16 

1.17(1) 

- 


- 



32 

1.40(-2) 

3.06 

8.14(-3) 

3.31 


MP5 

64 

5.05(-4) 

4.80 

3.01C-4) 

4.76 

11.86 


128 

1.63(-5) 

4.96 

9.74(-6) 

4.95 



256 

5.25(-7) 

4.95 

3.14 (- 7) 

4.96 



16 

1-17(1) 

- 


- 



32 

1.40(-2) 

3.06 

8.14(-3) 

3.30 


Unlim. 

64 

5.05(-4) 

4.80 

3.01(-4) 

4.76 

6.18 


128 

1.63(-5) 

4.96 

9.71(6) 

4.95 



256 

5.25(-7) 

4.95 

3.14 (- 7) 

4.96 



Table 2(b): At/Ax — 0.4, CPU time for constant reconstruction: 1.05 sec. 


Scheme 

'V 

Loo error 

Loo> order 

L\ error 

L\ order 

CPU time - sec. 


16 

2.39(-l) 

- 

1.07 (-1) 

- 



32 

3.7 4(-2) 

2.68 

1.8 7 (- 2) 

2.52 


WEN 05 

64 

3.26(-3) 

3.52 

1.79(-3) 

3.39 

3.30 


128 

3.00(-4) 

3.44 

1 -11 (■ 4) 

4.01 



256 

1.25(-5) 

4.58 

6.17(-6) 

4.17 



16 

1-21(1) 

- 

8.01(-2) 

- 



32 

1.7 7 (- 2) 

2.77 

1.03(-2) 

2.96 


MP5 

64 

1.10(-3) 

4.01 

G. 15 (- 4) 

4.06 

2.02 


128 

9.50(-5) 

3.54 

5.05(-5) 

3.61 



256 

1.04(-5) 

3.19 

5.42(-6) 

3.22 



16 

1.21(-1) 

- 

OO 

o 

t— 1 

• i 

to i 

- 



32 

1.7 7 (- 2) 

2.77 

1.03(-2) 

2.96 


Unlim. 

64 

1.10(-3) 

4.01 

6.17(-4) 

4.06 

1.29 


128 

9.50(-5) 

3.54 

5.04(-5) 

3.61 



256 

1.04(-5) 

3.19 

5.42(-6) 

3.22 
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4.2 Advection of a piecewise continuous function 

Next, the initial condition is given by 

n 0 (x) = exp( —log(2)(x + 0.7) 2 /0.0009) if - 0.8 < x < -0.6, 

u 0 (x) =1 if - 0.4 < x < -0.2, 

uo(ff) = 1 — |10(a: — 0.1)| if 0 < x < 0.2, (4-1) 

u 0 (x) = [1 - 100(3- - 0.5) 2 ] 1/2 if 0.4 < 3 - < 0.6, 

^ 0 ( 3 :) = 0 otherwise. 

This initial condition includes a Gaussian wave, a square wave, a triangular wave, and a semi¬ 
ellipse. We use 200 cells with a = 0.4. The solutions at t = 2 (after one period or 200 cells) and 
t = 20 (ten periods) are shown in Figs. 4.1 and 4.2 respectively. The solid line represents the exact 
solution. Also shown are the computing times of the various schemes. Again, note that the MP5 
solutions compare well with those by EN03 and WEN05 schemes. 

Resolution at discontinuities can be enhanced by using steepening techniques as in [4], [21]. and 
[6]. These techniques are expensive and, while they are effective in one dimension, it is still not 
clear how well they perform in multi-dimensions. Here, we will limit our study to the base schemes 
only. 

4.3 Euler system in one dimension 

In the following three problems, the CFL number is 0.4, and the spatial domain is [-1.1]. 
For the initial conditions, unless otherwise stated, the subscript /, denotes -1 < x < 0. and R. 
0 < x < 1. The final time is denoted by tj, and the total number of cells, by N. 

1. Sod’s problem [13] 

(Pl-ul,Pl) =( 1 , 0 , 1 ), 

[piu UR,p R ) = (0.125,0,0.1). 
tf = 0.4, 

N = 100. 

Since this problem starts from a singularity, smaller time steps are used initially as described in 
[6]. The density field from the MP5 scheme is shown in Fig. 4.3. Note that the contact discontinuity 
and the shock are resolved with high resolution. 
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Figure 4.1: Advection over a period bv (a) MP5, (b) WKN05, (c) EN03, (d) Unlimited sclieme 
A//Ax = 0.4, Ax = 2/200, t — 2, CPU time for constant reconstruction = 0.63 sec. 
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Figure 4.2: Advection over 10 periods by (a) MP5, (b) WEN05, (c) EN03, (d) Unlimited scheme 
At/Ax = 0.4, Ax = 2/200, t = 20, CPU lime for constant, reconstruction = 2.69 sec. 
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2. Lax’s problem [10] 


(pl,u l , Pl ) = (0.445,0.698,3.528), 

(pr,ur,pr) = (0.5,0,0.571), 
t f =0.32, 

N = 100. 

The density field from the MP5 scheme is shown in Fig. 4.4. Again, the contact discontinuity 
and the shock are well-resolved. 

3. Shu’s problem [14] 

In this problem, a moving shock wave interacts with a density disturbance and generates a flow 
field with both smooth structure and discontinuities. Here, L stands for — 1 < x < —0.8, and R , 
— 0.8 < x < 1 . The initial conditions, final time, and number of mesh points are 

(pl*ul,Pl) = (3.857143,2.629369,10.3333), 

(PR, ur, pr) = (1 + 0.2 sin( 57 T 2 ’), 0, 1), 
t f =0.36, 

N = 300. 

Since the exact solution is not known, the solution by EN03 with 800 cells is used in its place. 

The results of MP5 and WEN05 are shown in Fig. 4.5. The MP5 scheme captures the shock 
with high resolution and resolves all local extrema accurately. 

4.1 Euler system in two dimensions 

We present results for two well known problems. 

1 . Oblique shock reflection [1] 

The domain [0,4] x [0, 1] is covered by a uniform mesh of 60 x 20 cells. The boundary conditions 
are: at the bottom, solid boundary; at the right, supersonic outflow; at the left, the conditions are 
fixed with 

(p,u,v,p) = ( 1 , 2 . 9 , 0 , 1 / 7 ); 

and at the top, 

(p, u, v,p) = (1.69997,2.61934, -0.50632,1.52819). 

Under these conditions, an oblique shock forms from the top left corner and is reflected by the 
bottom boundary. Initially, flow conditions at the left boundary are set throughout the whole 
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domain. After 10,000 iterations, the solution essentially reaches steady state. The residua] drops 
roughly 2 orders of magnitude for MP5, WEN05, and EN03, while for the minmod and first-order 
upwind schemes the residual drops to machine zero. (Note that the minmod scheme is defined by 
w i /2 = w o + |minmod (w, - w 0 , w 0 - w^), where w is the characteristic flux in this case.) 

The pressure along the line y = 0.55 (j=ll) is shown in f ig. 4.6. Concerning accuracy, it 
can be seen that the higher-order schemes have small oscillations about the exact solution. These 
oscillations are reduced on finer grids for all three schemes. Note that the MP5 scheme yields a 
highly accurate solution. 

The computing times of the various schemes are also shown in Fig. 4.6. The first-order and 
minmod schemes are coded here with the local characteristic decompositions over the full five-point 
stencil. This first-order scheme represents the most efficient reconstruction in this framework, and 
its CPU time reflects the overhead of the characteristic decomposition. It can be seen from Fig. 4.6 
that the overhead is more than half of the total computing time. In other words, unlike the case of 
advection, the computing time of the reconstruction step for the Euler equations is less than one 
third of the total time. 

2. Double. Mach reflection [20] 

The computational domain is [0, 4] x ( 0 , 1]. The reflecting wall is from (1 /6,0) to (4.0). Initially, 
a Mach 10 shock is incident on this wall at (1/6,0) making an angle of sixty degrees with the .r-axis. 
To the right of the shock is undisturbed fluid of uniform pressure 1 and density 1.4. To the left of 
the shock, the conditions are 

(p.u.ihp) = (8.0,7.1447,-4.125,116.5). 

As the shock reflects off the wall, a diffraction pattern is formed. The final time is tj = 0.2. A 
detailed description of the problem and various solutions can be found in [20], 

The boundary conditions are: at the bottom, from (0,0) to (1/6,0), linear extrapolation; from 
(1/6,0) to (4,0), solid boundary; at the right, linear extrapolation; at the left, supersonic inflow; 
at the top, time-dependent conditions determined by the exact motion of the Mach 10 shock. 

The MP5 and WEN05 solutions, obtained using a 210 x 60 mesh, are shown in Fig. 4.7. It can 
be seen that both schemes capture all the significant features of the solution such as the two Mach 
stems and the wall jet. 
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5. Conclusions 


A new class of high-order schemes for the numerical solution of hyperbolic conservation laws was 
introduced. The key feature of these schemes is the reconstruction procedure which combines an 
accurate interface formula with a monotonicity-preserving constraint. The constraint is designed 
to preserve accuracy and to work well with Runge-Kutta time stepping. It is shown that, for 
advection, if the data are monotone, then the solution is also monotone under a time step restriction. 
Numerical experiments confirm that the resulting scheme is accurate in smooth regions, resolves 
discontinuities with high resolution, and is also efficient. The new scheme compares favorably with 
state-of-the-art schemes such as EN03 and WEN05. 

Appendix-A 


In this appendix, we give the listings of the three higher-order reconstruction algorithms in For¬ 
tran. V(J) are the cell averages vj and VL(J) are the computed interface values ^ +1 / 2 * DMM(X,Y) 
is the minmod function of two arguments while DM4(W,X,Y,Z) is the minmod function of four 
arguments. 


c-MP5 RECONSTRUCTION 

c-DMM(X.Y) = 0.5*( SIGN( 1. ,X J + SICN( 1 , Y ))* 

k M1N ( A B S ( X ), A B S ( Y )) 

c D M 4 (W, X, Y, Z) - 0.125*( 51G N (1., W ) + SIGN(1..X))* 

k ABS( ( SIGN( 1., W ) + SIGN{ 1 -,Y ))* 

k ( SIGN(1., W ) + SIGN( 1. ,Z )) ) 

At *M IN( A BS( W),ABS(X ), A B S( Y ), A BS(Z)) 

c 

B 1 = 0.0166(56666667 
B2 = 1 333333333333 

ALPHA = 4 
EPS M = IE-10 

c 

VOR = Bl*(2.*V(J-2)-13.*V(J-l) 
k -j- 47.*V(J) + 27.* V( J -j-1) 
k - 3. 1, ‘V(J-j-2) ) 

VMP= V(J) -f DMM(V(J+1)-V( J),ALPHA*(V( J)-V( J-l))) 
IF((VOR*V(J))*(VOR-VMP).LE EPSM) THEN 
VL( J ) = VOR 
ELSE 
ELSE 
C 

DJM1 = V( J-2)-2 *V( J-l)+ V(J ) 

DJ = V( J-1)-2.*V(J )+ V(J + 1) 

DJP1 = V(J )-2.*V( J+l)4 V(J + 2) 

DM4JPH= DM4( 4 *DJ-DJP1,4 "D JP 1-D J,D J ,D JP 1) 
DM4JMH= DM4(4.*DJ-DJM1,4.*DJM1-DJ,DJ 1 DJM1) 

VUL = V(J) + ALPHA *(V(J)-V( J-l)) 

VAV = 0.5*(V(J) + V(J + 1)) 

VMD = VAV - 0.5*DM4JPH 

VLC = V(J ) -f 0.5*(V( J)-V(J-l)) + B2*DM4 JMH 
VMIN = MAX(MIN(V(J) f V(J + l’),VMD), 
k MIN ( V (J ), V U L, V L C)) 

VM AX = MLN(MAX(V(J),V(J + 1),VMD), 
k MAX(V(.I),VUL,VLC)) 

VL( J ) = VOR 4 DMM( VMIN-VOR.VMAX-VOR) 

END1F 

C-WEN05 RECONSTRUCTION 


21 



EPS W = l.E-6 

B1 = 1.083333333333 

B2 = 0.166666666667 

DJM1 = V( J-2)-2 *V( J-l) + V(J ) 

EJMl = V(J-2)-4 *V(J-l) + 3 *V(J ) 

DJ = V( J-l )-2.*V( J ) + V(J+1) 

EJ = V(J-l)- V(J+1) 

DJP1 = V(J )-2 *V(J+1) + V(J-f2) 

EJP1 = 3*V( J)-4 *V(J + 1)+V( J + 2) 
c 

DISO = B1*DJM1*DJM1 + 0.2S“EJM J*EJMl -f EPSW 
DIS1 = B1*DJ*DJ + 0.25*EJ*EJ + EPSW 
DIS2 = Bl*DJPl*DJPl + 0 25*EJPl*EJPl + EPSW 
C 

Q30 = 2 *V(J-2)-7.*V(J-l) + 11 *V(J ) 

Q31 = - V( J-l ) + 5.*V( J )+ 2.* V( J-f 1) 

Q32 = 2*V(J ) + 5.*V(J + l) -V(J + 2) 

C 

D01 = DIS0/DIS1 

D02 = DIS0/DIS2 

A1BA0 = 6*D01*D01 

A2BA0 = 3.*D02*D02 

WO = 1/(1 + A1 B AO + A2BA0) 

W1 = A1 B AO* WO 
W2 = 1 - WO - W1 

VL(J) = B 2*( W0*Q30 + W1*Q31 + W2*Q32 } 

C 

C-EN03 RECONSTRUCTION 

DATA CM(l,l),CM(l,2),CM(l,3)/2 ,-7 ,11./ 

DATA CM(2,l) 1 CM(2 1 2),CM(2 t 3)/-l. t V,2./ 

DATA CM(3,1),CM{3,2),CM( 3,3)/2 ,5 1 / 

B1 = 0 166666666667 

SP = A B5( V(J + 1) - V( J ) ) 

SM = AB$( V(J ) - V(J-l) ) 

DJ = A BS( V( J + l) - 2 *V( J > + V(J-l) ) 

C 

IF(2.“5P .GT. SM) THEN 

DJM1 = ABS( V(J )- 2.*V(.M) -f V(J-2) ) 
IF(DJGT2*DJMl) THEN 
ID = 1 
ELSE 
ID = 2 
ENDIF 
ELSE 

DJP1 = A BS( V( J + 2)-2.*V{ J + l) + V(J ) ) 

IF ( 2. * D J P1 GT. DJ) THEN 
ID = 2 
ELSE 
ENDIF 
ENDIF 
C 

VL(J)= { CM(ID,1 )*V( J-3+ID) + CM{ ID,2 )" V( J-2 + ID ) + 
& CM(ID,3)*V( J-l + ID) )*B 1 
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Figure 4.5: Shu's problem with MP5 and WEN05 schemes. 
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Figure 4.6: Oblique shock reflection problem with various schemes, a — 0.4, 60 X 40 grid 

CPU time in seconds is shown above each plot. 
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Figure 4.7(a): Double Mach reflection problem with MP5. CFL = 0.4, 240 X 60 grid. 
30 density contours from 1.73 to 21. CPU time = 5553 seconds. 

CPU time for constant reconst ruction = 4122 seconds. 
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Figure 4.7(b): Double Mach reflection problem with YVEN05, CFL = 0.4, 240 X 60 grid. 
30 density contours from 1.73 to 21. CPU time = 6843 seconds. 

CPU time for constant recons1ruction= 4122 seconds. 
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